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We demonstrate by explicit calculations in 3+1 dimensional viscous relativistic fluid dynamics 
how two-particle pseudorapidity correlation functions in heavy ion collisions at the LHC and RHIC 
depend on the number of particle producing sources and the transport properties of the produced 
medium. In particular, we present results for the Legendre coefficients of the two-particle pseu¬ 
dorapidity correlation function, in Pb+Pb collisions at 2760 GeV and Au+Au collisions at 

200 GeV from viscous hydrodynamics with three dimensionally fluctuating initial conditions. Our 
results suggest that the Un.m provide important constraints on initial state fluctuations and the 
transport properties of the quark gluon plasma. 


INTRODUCTION 

The study of the Fourier coefhcients Vn of the az¬ 
imuthal particle distribution in heavy ion collisions at 
the Relativistic Heavy Ion Collider (RHIC) and the Large 
Hadron Collider (LHC) has led to great insight into the 
initial state fluctuations and hydrodynamic evolution of 
the produced system [1]. In particular, it has led to the 
conclusion that the system has a very small viscosity, 
close to the lower bound conjectured using AdS/CFT 
correspondence [2]. 

Recently, the ATLAS collaboration has presented re¬ 
sults on the expansion of two-particle pseudorapidity cor¬ 
relations into Legendre polynomials [3]. The obtained 
coefficients contain important information on the fluctu¬ 
ations of the particle multiplicity in pseudorapidity. 

In this letter we explore how two-particle pseudorapid¬ 
ity correlations from hydrodynamic calculations can give 
insight into the number of sources for particle produc¬ 
tion and their correlations, as well as the shear and bulk 
viscosity of the system. To achieve this, we introduce a 
simple initial state model that extends the conventional 
Monte Carlo (MC) Glauber model [4] to include fluctu¬ 
ations of valence quarks in rapidity and thus produces 
three dimensional fluctuating distributions of net baryon 
number and entropy density. We then use this model to 
generate the input for 3+1 dimensional viscous hydrody¬ 
namic calculations and compute rapidity distributions of 
charged hadrons and two-particle rapidity correlations. 
We then analyze the effect of i) the number of sources 
and ii) the shear and bulk viscosity of the system on the 
Legendre coefflcients. 

INITIAL STATE MODEL AND 
HYDRODYNAMIC EVOLUTION 


Gated string with many parameters, we introduce a sim¬ 
ple model that is a straight forward extension to the MC 
Glauber model to include longitudinal fluctuations. In 
particular, we will employ a Monte Carlo implementation 
of the Lexus model [9] to provide a simple parametriza- 
tion of the rapidity distribution of wounded nucleons (or 
constituent quarks), which is based on experimental ob¬ 
servations in proton-proton collisions. 

3DMC-Glauber model For each nucleus we sample the 
three-dimensional spatial distribution of nucleons from 
a Woods-Saxon distribution. We then sample valence 
quark positions within each nucleon from an exponential 
distribution. After overlaying the two nuclei in the trans¬ 
verse plane, separated by the sampled impact parameter 
b, wounded quarks are determined using the quark-quark 
inelastic cross-section aqq, which can be obtained geomet¬ 
rically by requiring that the experimentally determined 
nucleon-nucleon cross section is recovered. We do, how¬ 
ever, treat the quark-quark cross section in heavy ion 
collisions as an independent, free parameter. When we 
choose aqq = 3 mb at ^/s = 200 GeV, we can reproduce 
the experimental multiplicity and pseudorapidity distri¬ 
bution of charged hadrons without including additional 
negative binomial fluctuations. We employ a Gaussian 
wounding where quarks are determined to be participants 
with a Gaussian probability [10, 11]. 

As mentioned above, to obtain the longitudinal distri¬ 
bution of the initial baryon number we employ the Lexus 
model [9]. Here, the idea is that the rapidity distribution 
of nucleons in heavy ion collisions can be obtained by lin¬ 
ear extrapolation from the distribution in p+p collisions. 
That distribution is parametrized and fit to experimental 
data. We use this model for valence quarks such that the 
probability for a quark with rapidity yp to end up with 
rapidity y after a collision with a quark with rapidity yp 
(from the other nucleus) is given by [9]: 


Fluctuations in rapidity have been included in hydro- 
dynamic calculations via UrQMD [5, 6], EPOS [7], or 
AMPT [8] initial conditions. To study the effect of fluc¬ 
tuations in rapidity in addition to fluctuations in the 
transverse plane without being dependent on a compli- 


Q{y - yT,yp -yT,y- yp) 
, cosh(y - yp) 
smh{yp - yp) 


+ (1 - X)S{y - yp ), 


( 1 ) 


which corresponds to a flat distribution in longitudinal 
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momentum ■ In the original work [9] A is the fraction 
of nucleon-nucleon scatterings that result in a hard col¬ 
lision. Here, we treat A as a free parameter, regulating 
the stopping power of the collision. It can be fixed by 
fitting the net baryon distribution to experimental data. 
Generally, we find a good fit to the net baryon rapidity 
distribution when aqqX ~ 2 mb. 

The initial rapidities are distributed according to the 
quarks’ x value, which is determined by the valence quark 
parton distribution functions (PDFs).^ Initial rapidities 
are thus y = ±ln(a;-ys/2mg), with the sign depending on 
whether it is a quark in the projectile (-I-) or the target 
(-). We employ a valence quark mass of = 0.34 GeV 
and assume zero transverse momentum initially. 

To systematically organize the collisions of all quarks, 
we have quark pairs collide subsequently with increas¬ 
ing inter-quark distance I\z. We then convert rapidity 
to space-time rapidity ps to obtain a three dimensional 
event-by-event distribution of quarks. To assign a baryon 
density, this distribution needs to be smeared and we 
do this by introducing Gaussians in the transverse plane 
with width ctt = 0.4 fm, and width in space-time rapidity 
of (T = 0.2. 

Next we determine the entropy density distribution. 
We deposit a fixed entropy between the rapidities of each 
wounded quark and its last collision partner and smear 
the edges with half Gaussians of the same width as used 
for the baryon density distribution. This leads to the 
following form of the rapidity dependence of the entropy 
density per wounded quark pair 

s{y,yp, yr) = Afexp[-e{-y + mm{yp, yr)) 

X {y- min(yp, 2 /t))^/2(t^ 

- e{y - max{yp,yT)) 

X {y- max( 2 /p, yT)f/ 2 a'^] 

X (max{yp, yp) - min(yp, yp) + , (2) 

where yp and yp are the rapidities of the colliding quarks, 
and Af determines the normalization, which is the same 
for every “tube” and adjusted to fit the experimental 
data. 

In the transverse plane, we smear the entropy den¬ 
sity around the center of mass position of each pair by a 
Gaussian of width ap = 0.4 fm. 

This way of initializing the entropy density leads to 
the correct centrality dependence of the total multiplic¬ 
ity. We note that energy and momentum conservation is 
not explicitly fulfilled, however, we are not including any 
transverse momentum production or very high momen¬ 
tum partons, which will carry away some of the energy 


^ The X values are sampled from CTIO NNLO parton distribution 
functions [12] at = 1 GeV^ with EPS09 nuclear correction 
[13] using LHAPDF 6.1.4 [14]. 


and momentum of the collision and are not part of the 
bulk medium we are describing. 

A different initial state model using random rapidities 
for wounded nucleons was employed in [15]. 

Hydrodynamics We integrate the above initial condi¬ 
tion into the 3-1-1 dimensional viscous relativistic fluid 
dynamic simulation Music [16-19]. In addition to nu¬ 
merically solving the equations for the conservation of 
energy and momentum = 0, and the baryon cur¬ 

rent dfiJp = 0, we solve the relaxation-type equations 
derived from kinetic theory [20, 21] 

unH-b n = —— i5nnII0-I-(3) 

T^Ti{yix) + = 2pa^'^ - 

- . (4) 

The transport coefficients Tn, (5nn, Anu-j t-k, <5^71-, <^7^ T'-n-K, 
and A.^-!! are fixed using formulas derived from the Boltz¬ 
mann equation near the conformal limit [21]. At zero 
baryon chemical potential the ratio p/s will be chosen 
to be constant in this work, and the temperature depen¬ 
dence of the ratio of bulk viscosity to entropy density 
C/s is parametrized as in [22, 23], except that we grad¬ 
ually reduce the constant value at low temperature such 
that effects of the bulk Sf corrections [24] are kept to 
a minimum. Because we include finite baryon chemical 
potential yp > 0, we replace s in the above expressions 
by (e + P)/T, motivated by the relativistic limit of the 
fluidity measure introduced in [25]. 

The equation of state, which includes finite baryon 
chemical potential, is constructed by interpolating the 
pressures of hadronic resonance gas and lattice QCD 
[26, 27]. 

We leave a detailed description of the initial state 
model and the newly constructed equation of state for 
a longer paper in the future. 

After Cooper-Frye freeze out at an energy density of 
O.lGeV/fm^, the calculation of thermal spectra includ¬ 
ing Sf corrections [24] for shear and bulk viscosities, 
and resonance decays, we obtain the final hadron spec¬ 
tra as functions of transverse momentum pp and pseudo¬ 
rapidity Pp. 


TWO PARTICLE RAPIDITY CORRELATIONS 


The Pp integrated {pp > 0.5 GeV) event-by-event ra¬ 
pidity distributions dN/dpp are then used to determine 
the two-particle rapidity correlation function [28] 


C{pi,P2) 


{N{p,)N{p2)) 

{N{p,)){N{p2)) ’ 


( 5 ) 


where N{p) = dN/dpp. We follow the ATLAS collab¬ 
oration [3] and remove the effect of residual centrality 
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dependence in the average shape {N{r])) by redefining 
the correlation function as [29] 


CNijihm) 


Cp{m)Cp{ri2) ’ 


( 6 ) 


reduction of the longitudinal pressure keeps the pseudo¬ 
rapidity distribution more compact as shown previously 
[18, 37], leading to smaller suppression of the lower 
that are sensitive to long range correlations. 


where 


1 f 

Cy{rixi2) = — J ^C{r]i,r]2)dr]2/i- (7) 

Importantly, the resulting distribution is then normalized 
such that the average value of C'Ar(77i, 772 ) is one. 

Following [29, 30] C'Ar( 77 i, 772 ) is then decomposed into 
Legendre polynomials. The Legendre coefficients are 
given by 

^n,m ~ J 

T'„(77i)Tm(772) -I- Tn{ri2)Tm{m) dm dm 


Pb+Pb 2760 GeV 
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initial entropy density 
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•; 
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(1,1)(2,2)(3,3)(4,4)(5,5)(6,6)(1,3)(2,4)(3,5)(4,6)(5.7) 


where Tn{m) = + -P" stan¬ 

dard Legendre polynomials. The an,m are related to the 
Legendre coefficients of the single particle distribution 
a„ via anm = {o-nCLm), where the an are defined through 

Niv)/{Nm=l+En<^nTn{v) [3]- 

RESULTS 

For the calculations presented in this work, the pa¬ 
rameters of the initial state model are adjusted to fit 
the measured dN/drjp [31-33] and dN/d'^px [34] distri¬ 
butions of charged hadrons, as well as the net-baryon ra¬ 
pidity distribution [35, 36], when using the T-dependent 
bulk viscosity and 77 /s = 0.12. In particular we use a 
quark-quark cross section of 3 mb and A = 0.66. The 
hydrodynamic evolution starts at tq = 0.38 fm/c. 

We present results for the of the fluctuat¬ 

ing initial entropy density distribution (averaged over 
the transverse directions) and the final obtained 

from the produced charged hadrons with |r 7 i|, I 772 I < 2.4 
for 20-25% central Pb-|-Pb collisions at 2760 GeV in 
Fig. 1. Apart from the = (a^), the next largest cor¬ 
relations are those for m = n + 2. For this combination 
of Legendre indices, the an,m are negative, in qualitative 
agreement with the experimental data from ATLAS [3]. 
The an,n+i vanish in symmetric collision systems. 

The final y^|an,m| show a reduction relative to those 
obtained from the initial entropy density distribution. 
This reduction increases with 71 , 771 , showing that the hy¬ 
drodynamic expansion smears out the shorter range cor¬ 
relations more efficiently. The effect of viscosity is to 
slightly increase the lower an,m and decrease the higher 
ones. This can be understood as the effect of diffusion 
that destroys shorter scale structures in N{rip), while the 


FIG. 1. (Color online) Legendre coefficients labeled 

by {n, m) from the initial entropy density distribution and the 
final charged hadrons for 20-25% central collisions. Hydrody¬ 
namic expansion reduces the coefficients more efficiently with 
increasing n, m. Including shear and bulk viscosity leads to 
an increase of the lower \J\an,m.\ and a decrease of the higher 
\/|nn,m I. 

Generally, the n dependence of the calculated final 
state ^/\an,m\ is stronger than observed in the experi¬ 
mental data [3]. The calculated an,m agree with the ex¬ 
perimental data for small n = 777 = 1, but are smaller 
than the experimental values for all other 77 , 777 . 

The initial state values of \J\an,m\ agree fairly well 
with the experimental data, which is interesting, how¬ 
ever, a direct comparison is difficult in this case. 

The final calculated for large 77,777 are much 

smaller than the data for any choice of transport param¬ 
eters. Since we are only interested in the effects of flow 
in this study we have not included any non-flow effects, 
such as short range correlations from resonance decays 
(see e.g. [38]). The large difference between the calcula¬ 
tion and the experimental data for larger 77 , m indicates 
that these have a significant effect for 77,777 > 3. We 
leave a quantitative study of this effect using statistical 
hadronization for future work. 

The same effect is seen in 65-70% central collisions 
shown in Fig. 2. Here the coefficients extracted from the 
initial distribution already underestimate the experimen¬ 
tal data and after hydrodynamic evolution most y/\an,m\ 
are largely underestimated. 

The agreement of the lower an,m with experimental 
data for 20-25% centrality, but disagreement for 65-70%, 
indicates that the simple initial state model does not de¬ 
scribe correctly the centrality dependence of the number 
of sources, which dictates the degree of fluctuations. 
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FIG. 2. (Color online) Legendre coefficients y/\an,m.\ from 
the initial entropy density distribution and the final charged 
hadron distribution for 65-70% central collisions. 


In Fig. 3 we show predictions for y^\an,m\ in Au-|-Au 
collisions at 200 GeV. We find an overall increase of all 
compared to the Pb-bPb result at 2760 GeV by 
approximately 50 %. We further study the effect of vary¬ 
ing the shear viscosity to entropy density ratio ry/s from 
0.12 to 0.2, and the effect of bulk viscosity separately. As 
opposed to the change of ry/s from 0 to 0.12 for Pb-|-Pb 
collisions (see Fig. 1), we generally find an increase (or 
no change) of all y/\an,m\ when increasing rj/s. This in¬ 
dicates that the reduction of the longitudinal pressure 
dominates over the increased diffusion. Including bulk 
viscosity has a similar effect as increasing the shear vis¬ 
cosity, however, the effect is not large compared to the 
statistical errors of our results. 
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FIG. 3. (Golor online) Legendre coefficients \/\an,m.\ from 
the initial entropy density distribution and the final charged 
hadron distribution for 20-25% central Au-|-Au collisions at 
200 GeV. Larger shear viscosity leads to larger ^^/\an,m\■ So 
does the inclusion of bulk viscosity. 


Next we explicitly study the dependence of the an,m 
on the number of sources in the initial state. Fig. 4 shows 
how the of the initial entropy density distribution 
are increased when reducing the number of sources by 
employing nucleons instead of valence quarks. As alluded 
to above, this effect is expected, because fewer sources 
lead to more fluctuations. 
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FIG. 4. (Golor online) Legendre coefficients from the 

initial entropy density distribution. More sources lead to less 
fluctuations of the entropy density in pseudorapidity and thus 
smaller a^.m. 


Finally, in Fig. 5 we present a prediction for ^J\an,m\ 
extracted from the net-baryon distribution. We see that 
the hydrodynamic evolution does not dampen the higher 
an,m more than the lower ones as was the case for the 
charged hadron distribution. We have checked that the 
ideal results agree with the ones shown for rj/s = 0.12 
and {C/s){T) within errors. We argue that measuring 
these quantities e.g. at RHIG and comparing to various 
models could reveal important details on the mechanism 
of baryon stopping and transport in heavy ion collisions. 


CONCLUSIONS AND OUTLOOK 

We have presented the first calculation of two-particle 
pseudorapidity correlations in a viscous fluid dynamic 
framework and studied the dependence of our results 
on the initial state and the transport properties of the 
medium. We have shown that the number of particle pro¬ 
ducing sources affects the Legendre coefficients a„,m, as 
does the presence of hydrodynamic evolution. The shear 
and bulk viscosities of the medium also modify the 
coefficients. The effect of viscosity is two-fold. On the 
one hand, it reduces the longitudinal pressure, leading to 
an increase of the an^m- On the other hand, the diffu¬ 
sion has the opposite effect. This effect becomes larger 
for shorter range correlations, i.e., larger n,m. The un¬ 
derestimation of the experimental y^\an,m\ coefficients in 
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FIG. 5. (Color online) Legendre coefficients y/\an,m\ of the 
net-baryon distriution in 20-25% central 200 GeV Au+Au 
collisions at RHIC. In contrast to charged hadron correla¬ 
tions,the second coefficient -^|a( 2 , 2 ) | = \/ {a%) is larger than 
the first y'|a(i_i)|. 


particular for large n and m indicates that short range 
non-flow correlations have a large effect on these observ¬ 
ables. 

Based on our results we conclude that the study of 
Qn.m together with the Fourier coefficients Vn at various 
collision energies has the strong potential to constrain 
the initial state particle production mechanism as well 
as the transport properties of the medium. 

We further propose to perform the same experimental 
analysis using net-baryon or net-proton distributions to 
constrain in detail the mechanism of baryon stopping and 
transport in heavy ion collisions. 

It will also be very interesting to analyze the two- 
particle pseudorapidity correlations using a more so¬ 
phisticated initial state model, such as the IP-Glasma 
[19, 39, 40] in the future. This will require an extension 
of the boost-invariant Glasma picture to three dimen¬ 
sions. 

Note: An independent study of the two-particle pseu¬ 
dorapidity correlations within a similar framework and 
leading to similar conclusions has appeared simultane¬ 
ously with our work [41]. 
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